
#include <windows.h>
#include <time.h>
#include <malloc.h>
#include <stddef.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "vol_lip.h"

#ifndef KAR_DEPTH
# define KAR_DEPTH      20
#endif

#ifndef KAR_MUL_CROV
# define KAR_MUL_CROV   30
#endif

#ifndef KAR_SQU_CROV
# define KAR_SQU_CROV   30
#endif

#if (!ILLEGAL)


/*Function prototypes for static internal functions.*/

static void zsubmul(
	long r,
	verylong a,
	verylong b
	);

static void zalphasubmul(
	long r,
	verylong a,
	verylong b
	);

#define zhalt(msg)

static void zaddls(
	verylong a,
	verylong b,
	verylong c
	);

static void kar_mul(
	verylong a,
	verylong b,
	verylong *c,
	long shi
	);

static void kar_sq(
	verylong a,
	verylong *c,
	long shi
	);

static long zxxeucl(
	verylong ain,
	verylong nin,
	verylong *invv,
	verylong *uu
	);

#ifdef DOUBLES_LOW_HIGH
#define LO_WD 0
#define HI_WD 1
#else
#define LO_WD 1
#define HI_WD 0
#endif

#ifdef SINGLE_MUL

#define ExtractHiLo(Hi,Lo,x) \
{ \
double t=x+(4503599627370496.0); \
unsigned int *it = (unsigned int *)&t; \
Lo = it[LO_WD]; \
Hi = ((it[HI_WD]<<6)|(it[LO_WD]>>26)); \
Lo &= 0x3FFFFFF; \
Hi &= 0x3FFFFFF; \
}

#define zaddmulp(a,b,d,t) \
{ \
   double __lx = ((double) ((*a) + (*t))) + ((double) b)*((double) d); \
   register long  __lhi = 0, __llo = 0;\
   ExtractHiLo(__lhi,__llo,__lx);\
   (*a) = __llo;\
   (*t) = __lhi;\
}

#define zaddmulpsq(a,b,t) \
{ \
   double __lx = ((double) (*a)) + ((double) b)*((double) b); \
   register long  __lhi = 0, __llo = 0;\
   ExtractHiLo(__lhi,__llo,__lx);\
   (*a) = __llo;\
   (*t) = __lhi;\
}


typedef union { double d; unsigned long rep[2]; } d_or_rep;

#define FetchHiLo(hi,lo,x) \
{ \
   d_or_rep _l_xx; \
   _l_xx.d = x; \
   hi = _l_xx.rep[HI_WD]; \
   lo = _l_xx.rep[LO_WD]; \
}

static
void zaddmul(
        long ams,
        long *ama,
        long *amb
        ) 
{ 
        register long carry = 0; 
        register long i = (*amb++);
        register double dams = (double) ams;
        register double xx;
        register double yy;
        register unsigned long lo_wd, lo;
        register unsigned long hi_wd, hi;

        xx  =  ((double) (*amb))*dams + 4503599627370496.0;
        for (; i > 1; i--) 
        { 
                yy =  ((double) (*(++amb)))*dams +4503599627370496.0;
                FetchHiLo(hi_wd, lo_wd, xx);
                lo = lo_wd & 0x3FFFFFF;
                hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
                lo = lo + (*ama) + carry;
                *ama = lo & 0x3FFFFFF;
                carry = hi + (lo >> 26);
                ama++; 
                xx = yy;
        } 

        FetchHiLo(hi_wd, lo_wd, xx);
        lo = lo_wd & 0x3FFFFFF;
        hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
        lo = lo + (*ama) + carry;
        *ama = lo & 0x3FFFFFF;
        carry = hi + (lo >> 26);
        ama++; 
        *ama += carry; 
}


static
void zaddmulsq(
	long ams,
	long *ama,
	long *amb
	) 
{ 
        register long carry = 0; 
        register long i = ams;
        register double dams = (double) (*amb);
        double xx;
	register double yy;
        register unsigned long lo_wd,lo;
	register unsigned long hi_wd,hi;

	if (!ams) return;
        xx =  ((double) (*(++amb)))*dams + 4503599627370496.0;
        for (; i > 1; i--) 
        { 
                yy =  ((double) (*(++amb)))*dams + 4503599627370496.0;
		FetchHiLo(hi_wd, lo_wd, xx);
                lo = lo_wd  & 0x3FFFFFF;
                hi = ((hi_wd <<6)|(lo_wd >>26)) & 0x3FFFFFF;
                lo = lo + (*ama) + carry;
                *ama = lo & 0x3FFFFFF;
                carry = hi + (lo >> 26);
                ama++; 
                xx = yy;
        } 
	if (i==1)
	{
		FetchHiLo(hi_wd, lo_wd, xx);
        	lo = lo_wd  & 0x3FFFFFF;
        	hi = ((hi_wd <<6)|(lo_wd >>26)) & 0x3FFFFFF;
        	lo = lo + (*ama) + carry;
        	*ama = lo & 0x3FFFFFF;
        	carry = hi + (lo >> 26);
        	ama++; 
	}
        *ama += carry; 
}


static
void zsubmul(
	long ams,
	long *ama,
	long *amb
	) 
{ 
        register long carry = RADIX; 
        register long i = (*amb++);
        register double dams = (double) (RADIX-ams);
        double xx;
	register double yy;
        register unsigned long lo_wd,lo;
	register unsigned long hi_wd,hi;

        xx =  ((double) (*amb))*dams + 4503599627370496.0;
        for (; i > 1; i--) 
        { 
                yy =  ((double) (*(amb+1)))*dams + 4503599627370496.0;
		FetchHiLo(hi_wd, lo_wd, xx);
                lo = lo_wd & 0x3FFFFFF;
                hi = ((hi_wd<<6)|(lo_wd>>26)) & 0x3FFFFFF;
                lo = lo + (*ama) + carry;
                *ama = lo & 0x3FFFFFF;
                carry = hi + (lo >> 26);
                ama++; 
		carry += RADIXM - (*(amb++));
                xx = yy;
        } 
	FetchHiLo(hi_wd, lo_wd, xx);
        lo = lo_wd & 0x3FFFFFF;
        hi = ((hi_wd <<6)|(lo_wd >>26)) & 0x3FFFFFF;
        lo = lo + (*ama) + carry;
        *ama = lo & 0x3FFFFFF;
        carry = hi + (lo >> 26);
        ama++; 
	carry += RADIXM - (*amb);
        *ama += carry - RADIX; 
}

#elif KARAT

/*
	This definition of zaddmulp and zaddmulpsq assumes nothing
*/
#define zaddmulp(  a,  b,  d,  t) { \
        register unsigned long b1 = b & RADIXROOTM; \
        register unsigned long d1 = d & RADIXROOTM; \
        register unsigned long bd,b1d1,m,aa= (*a) + (*t); \
	register unsigned long ld = (d>>NBITSH); \
	register unsigned long lb = (b>>NBITSH); \
 \
        bd=lb*ld; \
        b1d1=b1*d1; \
        m=(lb+b1)*(ld+d1) - bd - b1d1; \
        aa += ( b1d1+ ((m&RADIXROOTM)<<NBITSH)); \
        (*t) = (aa >> NBITS) + bd + (m>>NBITSH); \
        (*a) = aa & RADIXM; \
}



#define zaddmulpsq(_a, _b, _t) \
{ \
	register long lb = (_b); \
	register long b1 = (_b) & RADIXROOTM; \
	register long aa = *(_a) + b1 * b1; \
 \
	b1 = (b1 * (lb >>= NBITSH) << 1) + (aa >> NBITSH); \
	aa = (aa & RADIXROOTM) + ((b1 & RADIXROOTM) << NBITSH); \
	*(_t) = lb * lb + (b1 >> NBITSH) + (aa >> NBITS); \
	*(_a) = (aa & RADIXM); \
}

#define zaddmul(ams, ama, amb) \
{ \
	register long lami; \
	register long lams = (ams); \
	register verylong lama = (ama); \
	register verylong lamb = (amb); \
	long lamcarry = 0; \
 \
	for (lami = (*lamb++); lami > 0; lami--) \
	{ \
		zaddmulp(lama, *lamb, lams, &lamcarry); \
	/* Be careful, the last lama is unnormalized */ \
		lama++; \
		lamb++; \
	} \
        *lama += lamcarry; \
}


#define zaddmulsq(sql, sqa, sqb) \
{ \
	register long lsqi = (sql); \
	register long lsqs = *(sqb); \
	register verylong lsqa = (sqa); \
	register verylong lsqb = (sqb); \
	long lsqcarry = 0; \
 \
	lsqb++; \
	for (; lsqi > 0; lsqi--) \
	{ \
		zaddmulp(lsqa, *lsqb, lsqs, &lsqcarry); \
		lsqa++; \
		lsqb++; \
	} \
	*lsqa += lsqcarry; \
/* Be careful, the last lama is unnormalized */ \
}

static void
zsubmul(
	long r,
	verylong a,
	verylong b
	)
{
	register long rd = RADIX - r;
	register long i;
	long carry = RADIX;

	for (i = (*b++); i > 0; i--)
	{
		zaddmulp(a, *b, rd, &carry);
		a++;
		carry += RADIXM - (*b++);
	}
	*a += carry - RADIX;	/* unnormalized */
}


#elif PLAIN

/*
	This definition of zaddmulp and zaddmulpsq assumes nothing
*/
#define zaddmulp(_a, _b, _d, _t) \
{ \
	register long lb = (_b); \
	register long ld = (_d); \
	register long b1 = (_b) & RADIXROOTM; \
	register long d1 = (_d) & RADIXROOTM; \
	register long aa = *(_a) + b1 * d1; \
 \
	b1 = b1 * (ld >>= NBITSH) + d1 * (lb >>= NBITSH) + (aa >> NBITSH); \
	aa = (aa & RADIXROOTM) + ((b1 & RADIXROOTM) << NBITSH) + *(_t); \
	*(_t) = ld * lb + (b1 >> NBITSH) + (aa >> NBITS); \
	*(_a) = (aa & RADIXM); \
}


#define zaddmulpsq(_a, _b, _t) \
{ \
	register long lb = (_b); \
	register long b1 = (_b) & RADIXROOTM; \
	register long aa = *(_a) + b1 * b1; \
 \
	b1 = (b1 * (lb >>= NBITSH) << 1) + (aa >> NBITSH); \
	aa = (aa & RADIXROOTM) + ((b1 & RADIXROOTM) << NBITSH); \
	*(_t) = lb * lb + (b1 >> NBITSH) + (aa >> NBITS); \
	*(_a) = (aa & RADIXM); \
}


#define zaddmul(ams, ama, amb) \
{ \
	register long lami; \
	register long lams = (ams); \
	register verylong lama = (ama); \
	register verylong lamb = (amb); \
	long lamcarry = 0; \
 \
	for (lami = (*lamb++); lami > 0; lami--) \
	{ \
		zaddmulp(lama, *lamb, lams, &lamcarry); \
	/* Be careful, the last lama is unnormalized */ \
		lama++; \
		lamb++; \
	} \
        *lama += lamcarry; \
}


#define zaddmulsq(sql, sqa, sqb) \
{ \
	register long lsqi = (sql); \
	register long lsqs = *(sqb); \
	register verylong lsqa = (sqa); \
	register verylong lsqb = (sqb); \
	long lsqcarry = 0; \
 \
	lsqb++; \
	for (; lsqi > 0; lsqi--) \
	{ \
		zaddmulp(lsqa, *lsqb, lsqs, &lsqcarry); \
		lsqa++; \
		lsqb++; \
	} \
	*lsqa += lsqcarry; \
/* Be careful, the last lama is unnormalized */ \
}

static void
zsubmul(
	long r,
	verylong a,
	verylong b
	)
{
	register long rd = RADIX - r;
	register long i;
	long carry = RADIX;

	for (i = (*b++); i > 0; i--)
	{
		zaddmulp(a, *b, rd, &carry);
		a++;
		carry += RADIXM - (*b++);
	}
	*a += carry - RADIX;	/* unnormalized */
}

#else

/*
	This definition of zaddmulp and zaddmulpsq presumes a two's complement
	machine in which integer overflow is ignored and where double
	precision arithmetic is fast. The 0.25 allows round to nearest
	or round towards zero (value being rounded should be integer
	except for roundoff.)

	DO NOT USE FOR alpha because of fradix_inv (unless NBITS is small)
*/

static double fradix_inv = 1.0 / RADIX;	/* Global constant */

#define zaddmulp(_a, _b, _d, _t) \
{ \
	register at = *(_a) + *(_t); \
	register long aa = (at + (_b) * (_d)) & RADIXM; \
 \
	*(_t) = (long) (0.25 + fradix_inv * (((double) (at - aa)) \
	               + ((double) (_b)) * ((double) (_d)))); \
	*(_a) = aa; \
}

#define zaddmulpsq(_a, _b, _t) \
{ \
	register at = *(_a); \
	register long aa = (at + (_b) * (_b)) & RADIXM; \
 \
	*(_t) = (long) (0.25 + fradix_inv * (((double) (at - aa)) \
				      + ((double) (_b)) * ((double) (_b)))); \
	*(_a) = aa; \
}


#define zaddmul(ams, ama, amb) \
{ \
	register long lami; \
	register long lams = (ams); \
	register verylong lama = (ama); \
	register verylong lamb = (amb); \
	long lamcarry = 0; \
 \
	for (lami = (*lamb++); lami > 0; lami--) \
	{ \
		zaddmulp(lama, *lamb, lams, &lamcarry); \
	/* Be careful, the last lama is unnormalized */ \
		lama++; \
		lamb++; \
	} \
        *lama += lamcarry; \
}


#define zaddmulsq(sql, sqa, sqb) \
{ \
	register long lsqi = (sql); \
	register long lsqs = *(sqb); \
	register verylong lsqa = (sqa); \
	register verylong lsqb = (sqb); \
	long lsqcarry = 0; \
 \
	lsqb++; \
	for (; lsqi > 0; lsqi--) \
	{ \
		zaddmulp(lsqa, *lsqb, lsqs, &lsqcarry); \
		lsqa++; \
		lsqb++; \
	} \
	*lsqa += lsqcarry; \
/* Be careful, the last lama is unnormalized */ \
}

static void
zsubmul(
	long r,
	verylong a,
	verylong b
	)
{
	register long rd = RADIX - r;
	register long i;
	long carry = RADIX;

	for (i = (*b++); i > 0; i--)
	{
		zaddmulp(a, *b, rd, &carry);
		a++;
		carry += RADIXM - (*b++);
	}
	*a += carry - RADIX;	/* unnormalized */
}
#endif


#ifdef ALPHA50
#define zaldivaddmulp(  a,  b,  d,  t) { \
        register long b1 = b & ALPHA50RADIXROOTM; \
        register long d1 = d & ALPHA50RADIXROOTM; \
        register long bd,b1d1,m,aa= (*a) + (*t); \
        register long ld = (d>>ALPHA50NBITSH); \
        register long lb = (b>>ALPHA50NBITSH); \
 \
        bd=lb*ld; \
        b1d1=b1*d1; \
        m=(lb+b1)*(ld+d1) - bd - b1d1; \
	aa += ( b1d1+ ((m&ALPHA50RADIXROOTM)<<ALPHA50NBITSH)); \
        (*t) = (aa >> ALPHA50NBITS) + bd + (m>>ALPHA50NBITSH); \
        (*a) = aa & ALPHA50RADIXM; \
}
#endif


#define zaddmulone(ama, amb) \
{ \
	register long lami; \
	register long lams = 0; \
	register verylong lama = (ama); \
	register verylong lamb = (amb); \
 \
	lams = 0; \
	for (lami = (*lamb++); lami > 0; lami--) \
	{ \
		lams += (*lama + *lamb++); \
		*lama++ = lams & RADIXM; \
		lams >>= NBITS; \
	} \
	*lama += lams; \
}

#define zalphaaddmulone(ama, amb) \
{ \
        register long lami; \
        register long lams = 0; \
        register verylong lama = (ama); \
        register verylong lamb = (amb); \
 \
        lams = 0; \
        for (lami = (*lamb++); lami > 0; lami--) \
        { \
                lams += (*lama + *lamb++); \
                *lama++ = lams & ALPHA50RADIXM; \
                lams >>= ALPHA50NBITS; \
        } \
        *lama += lams; \
}


/*
	zdiv21 returns quot, numhigh so

	quot = (numhigh*RADIX + numlow)/denom;
	numhigh  = (numhigh*RADIX + numlow)%denom;

Assumes 0 <= numhigh < denom < RADIX and 0 <= numlow < RADIX.
*/


#ifdef ALPHA

#define zdiv21(numhigh, numlow, denom, deninv, quot) \
{ \
	register long lq21 = (long) (((fradix * (double) (numhigh)) \
			 + (double) (numlow)) * (deninv)); \
	register long lr21; \
	register long dd; \
	long lprodhigh = 0; \
	long lprodlow = 0; \
 \
	zaddmulp(&lprodlow, lq21, denom, &lprodhigh); \
	dd = numhigh-lprodhigh; \
	lprodhigh=0; /* this line here because of compiler error */ \
	if (dd > 1) { \
		register long correction = (((double)dd)*(double)RADIX + numlow - lprodlow)*deninv; \
		lq21 += correction; \
		lprodlow = 0; \
		zaddmulp(&lprodlow, lq21, denom, &lprodhigh); \
		dd = numhigh-lprodhigh; \
	} \
	else if (dd < -1) { \
		register long correction = (((double)dd)*(double)RADIX + numlow - lprodlow)*deninv; \
		lq21 += correction; \
		lprodlow = 0; \
		zaddmulp(&lprodlow, lq21, denom, &lprodhigh); \
		dd = numhigh-lprodhigh; \
	} \
	lr21 = RADIX * dd + (numlow - lprodlow); \
	if (dd=lr21*(deninv)) { \
		lq21 += dd; \
		lr21 -= denom*dd; \
	} \
	if (lr21 < 0) \
	{ \
		while (lr21 < 0) { \
			lq21--; \
			lr21 += denom; \
		}; \
	} \
	else \
	{ \
		while (lr21 >= denom) \
		{ \
			lr21 -= denom; \
			lq21++; \
		}; \
	} \
	quot = lq21; \
	numhigh = lr21; \
}
#else

#ifdef ALPHA50

#define zdiv21(numhigh, numlow, denom, deninv, quot) \
{ \
        register long lr21; \
        register long lq21 = (long) (((alpha50fradix * (double) (numhigh)) \
                         + (double) (numlow)) * (deninv)); \
        long lprodhigh = 0; \
        long lprodlow = 0; \
 \
        zaldivaddmulp(&lprodlow, lq21, denom, &lprodhigh); \
        lr21 = ALPHA50RADIX *(numhigh - lprodhigh) + (numlow - lprodlow); \
        if (lr21 < 0) \
        { \
                do \
                { \
                        lq21--; \
                } while ((lr21 += denom) < 0); \
        } \
        else \
        { \
                while (lr21 >= denom) \
                { \
                        lr21 -= denom; \
                        lq21++; \
                }; \
        } \
        quot = lq21; \
        numhigh = lr21; \
}

#else

#ifndef WIN32
#define zdiv21(numhigh, numlow, denom, deninv, quot) \
{ \
        register long lr21; \
        register long lq21 = (long) (((fradix * (double) (numhigh)) \
                         + (double) (numlow)) * (deninv)); \
        long lprodhigh = 0; \
        long lprodlow = 0; \
 \
        zaddmulp(&lprodlow, lq21, denom, &lprodhigh); \
        lr21 = RADIX * (numhigh - lprodhigh) + (numlow - lprodlow); \
        if (lr21 < 0) \
        { \
                do \
                { \
                        lq21--; \
                } while ((lr21 += denom) < 0); \
        } \
        else \
        { \
                while (lr21 >= denom) \
                { \
                        lr21 -= denom; \
                        lq21++; \
                }; \
        } \
        quot = lq21; \
        numhigh = lr21; \
}
#else

#define zdiv21(numhigh, numlow, denom, deninv, quot) \
{ \
        register long lr21; \
        register long lq21 = (long) (((fradix * (double) (numhigh)) \
                         + (double) (numlow)) * (deninv)); \
 /* Following works in many two's complement architectures. */ \
        lr21 = (numhigh << NBITS) + numlow - lq21 * denom; \
        if (lr21 < 0) \
        { \
                do \
                { \
                        lq21--; \
                } while ((lr21 += denom) < 0); \
        } \
        else \
        { \
                while (lr21 >= denom) \
                { \
                        lr21 -= denom; \
                        lq21++; \
                }; \
        } \
        quot = lq21; \
        numhigh = lr21; \
}

#endif


#endif
#endif


#define MulLo(rres,a,b) \
{ \
   register double _x = ((double) a) * ((double) b); \
   double _t = _x+(4503599627370496.0); \
   register long *_lo = (long *)&_t; \
   rres =  _lo[LO_WD];\
}

extern void *calloc();

/* global variables */

/* for long division */
static double log16rad = -1.0;
static double epsilon = 0.0;
static double fradix = (double)RADIX;
static double fudge = -1.0;
#ifdef ALPHA
static double fudge2 = -1.0;
#endif
#ifdef ALPHA50
static double alpha50fudge = -1.0;
static double alpha50fradix = (double) ALPHA50RADIX;
#endif
/* for convenience */
static long oner[] = {1, 1, 1};
static verylong one = &oner[1];
/* for karatsuba */
static verylong kar_mem[5*KAR_DEPTH];
static long kar_mem_initialized = 0;

#ifndef WIN32
#include <sys/times.h>
#include <limits.h>
#ifndef CLK_TCK
#define CLK_TCK 60
#endif
#endif

void free_karmem ()
{
    register long aneg;

	if (kar_mem_initialized)
	{
		for (aneg = (5 * KAR_DEPTH) - 1; aneg >= 0; aneg--)
			zfree (&kar_mem[aneg]);
		kar_mem_initialized = 0;
	}
}

void
zsetlength(
	verylong *v,
	long len
	)
{
	verylong x = *v;

	if (x)
	{
		if (len <= x[-1])
			return;
#if (PRT_REALLOC>0)
			//fprintf(stderr,"%s reallocating to %ld\n", str, len);
			//fflush(stderr);
#endif
		x[-1] = len;
		if (!(x = (verylong)realloc((void*)(&(x[-1])), (size_t)(len + 2) * SIZEOFLONG)))
		{
			//fprintf(stderr,"%d bytes realloc failed\n", ((int)len + 2) * SIZEOFLONG);
			zhalt("reallocation failed in zsetlength");
		}
	}
	else if (len >= 0)
	{
		if (len < SIZE)
			len = SIZE;
#if (PRT_REALLOC>0)
			//fprintf(stderr,"%s allocating to %ld\n", str, len);
			//fflush(stderr);
#endif
		if (!(x = (verylong)calloc((size_t)(len + 2), (size_t)SIZEOFLONG)))
		{
			//fprintf(stderr,"%d bytes calloc failed\n", ((int)len + 2) * SIZEOFLONG);
			zhalt("allocation failed in zsetlength");
		}
		x[0] = len;
		x[1] = 1;
		x[2] = 0;
	}
	else
		zhalt("negative size allocation in zsetlength");
	x++;
	*v = x;
}

void
zfree(
	verylong *x
	)
{
	if (!(*x))
		return;
	{
		verylong y = (*x - 1);
		free((void*)y);
		*x = 0;
		return;
	}
}

void
zstart()
{
	double local_one = (double) 1;
	double local_half = 1 / (double) 2;
	epsilon = local_one;
	fudge = local_one + epsilon;

 /* the following loop is sick, but we have to: */
 /* comparing one to one + epsilon does not have */
 /* the desired effect on some machines, because */
 /* one + epsilon is Not cast to double (as it should) */
	while (local_one != fudge)
	{
		epsilon = epsilon * local_half;
		fudge = local_one + epsilon;
	}
	epsilon += epsilon;
#ifndef ALPHA_OR_ALPHA50
	if ((epsilon * RADIX) > local_one)
	{
		//fprintf(stderr,
		//	"zstart failure: recompile with smaller NBITS\n");
		//exit(0);
        return;
	}
#endif
	epsilon *= 3;
	fudge = fradix + epsilon * fradix;
#ifdef ALPHA
	fudge2 = fradix - epsilon * fradix;
#endif
#ifdef ALPHA50
	alpha50fudge = alpha50fradix + epsilon * alpha50fradix;
#endif
}


void
zzero(
	verylong *aa
	)
{
	if (!(*aa)) zsetlength(aa, (long)SIZE/*, "in zzero, first argument"*/);
	(*aa)[0] = (long) 1;
	(*aa)[1] = (long) 0;
}

void
zone(
	verylong *aa
	)
{
	if (!(*aa)) zsetlength(aa, (long)SIZE/*, "in zone, first argument"*/);
	(*aa)[0] = (long) 1;
	(*aa)[1] = (long) 1;
}

void
zcopy(
	verylong a,
	verylong *bb
	)
{
	register long i;
	verylong b = *bb;

	if (ALLOCATE && !a)
	{
		zzero(bb);
		return;
	}
	if (a != b)
	{
		if ((i = *a) < 0)
			i = (-i);
		zsetlength(&b, i/*, "in zcopy, second argument"*/);
		*bb = b;
		for (; i >= 0; i--)
			*b++ = *a++;
	}
}

void
zintoz(
	long d,
	verylong *aa
	)
{
	register long i = 0;
	register long anegative = 0;
	verylong a = *aa;

	if (!(*aa)) zsetlength(&a, (long)SIZE/*, "in zintoz, second argument"*/);
	*aa = a;
	if (d < 0)
	{
		anegative = 1;
		d = -d;
	}
	a[1] = 0;
	while (d > 0)
	{
		a[++i] = d & RADIXM;
		d >>= NBITS;
	}
	if (i > 0)
		a[0] = i;
	else
		a[0] = 1;
	if (anegative)
		a[0] = (-a[0]);
}

long 
ztoint(
	verylong a
	)
{
	register long d;
	register long sa;

	if (ALLOCATE && !a)
		return (0);
	if ((sa = *a) < 0)
		sa = -sa;
	d = *(a += sa);
	while (--sa)
	{
		d <<= NBITS;
		d += *(--a);
	}
	if ((*(--a)) < 0)
		return (-d);
	return (d);
}

char 
eulav(
	long i
	)
{
    if (i < 10)
        return (char)(i + '0');
    else if (i < 16)
        return (char)(i - 10 + 'A');
    else
        return '?';

/*	switch (i)
	{
		case 0:
			return ('0');
		case 1:
			return ('1');
		case 2:
			return ('2');
		case 3:
			return ('3');
		case 4:
			return ('4');
		case 5:
			return ('5');
		case 6:
			return ('6');
		case 7:
			return ('7');
		case 8:
			return ('8');
		case 9:
			return ('9');
		case 10:
			return ('A');
		case 11:
			return ('B');
		case 12:
			return ('C');
		case 13:
			return ('D');
		case 14:
			return ('E');
		case 15:
			return ('F');
		default:
			return ('?');
	} */
}

void zhswrite(
	char *str,
	verylong a
	)
{
	//static char *b = 0;
	//static bl = 0;
	char *b = 0;
	STATIC verylong aa = 0;
	register long i;
	register long cnt = 0;
	register long mb;
	register long lab = 0;

	if (!a)
		*str++ = '0';
	zcopy(a, &aa);
	if (aa[0] < 0)
	{
		aa[0] = -aa[0];
		*str++ = '-';
	}

	b = (char *)malloc((size_t)(aa[0] << 3));
	//if (!b)
	//	b = (char *)malloc((size_t)(bl = (aa[0] << 3)));
	//else if (bl < (aa[0] << 3))
	//	b = (char *)realloc((void*)b, (size_t)((bl = (aa[0] << 3)) * sizeof(char)));
	do
	{
		b[cnt] = eulav(aa[1] & 15);
		cnt++;
		zrshift(aa, (long) 4, &aa);
	} while ((aa[1] != 0) || (aa[0] != 1));
	if (!HEX_BLOCK) mb = 0;
	else mb = HEX_BLOCK - (cnt%HEX_BLOCK);
	if (mb == HEX_BLOCK) mb = 0;
	for (i = cnt; i--;)
	{
		*str++ = b[i];
		mb++;
		if (mb == HEX_BLOCK)
		{
			lab++;
			if (lab == HEX_BLOCKS_PER_LINE)
			{
				//if (i)
				//	fprintf(fn, "%c\n", OUT_LINE_BREAK);
				lab = 0;
				mb = 0;
			}
			else
			{
				mb = 0;
				//fprintf(fn, "%c", HEX_SEP_CHAR);
			}
		}
	}
	FREESPACE(aa);
    free (b);

    *str = '\0';
}

long 
zcompare(
	verylong a,
	verylong b
	)
{
	register long sa;
	register long sb;

	if (ALLOCATE && !a)
	{
		if (!b)
			return (0);
		if (b[0] < 0)
			return (1);
		if (b[0] > 1)
			return (-1);
		if (b[1])
			return (-1);
		return (0);
	}
	if (ALLOCATE && !b)
	{
		if (a[0] < 0)
			return (-1);
		if (a[0] > 1)
			return (1);
		if (a[1])
			return (1);
		return (0);
	}
 /* if (a==b) return (0); but who is comparing identical addresses? */
	if ((sa = *a) > (sb = *b))
		return (1);
	if (sa < sb)
		return (-1);
	if (sa < 0)
		sa = (-sa);
	a += sa;
	b += sa;
	for (; sa; sa--)
	{
		if (*a > *b)
		{
			if (sb < 0)
				return (-1);
			return (1);
		}
		if (*a-- < *b--)
		{
			if (sb < 0)
				return (1);
			return (-1);
		}
	/*
	 * depending on the relative speed of various operations it might be
	 * better to do the following: if ((aux=(*a--)-(*b--))>0) { if (sb<0)
	 * return (-1); return (1); } else if (aux<0) { if (sb<0) return (1);
	 * return (-1); } where aux is a register long
	 */
	}
	return (0);
}

void
znegate(
	verylong *aa
	)
{
	verylong a = *aa;

	if (!a)
		return;
	if (a[1] || a[0] != 1)
		a[0] = (-a[0]);
}

void
zsadd(
	verylong a,
	long d,
	verylong *b
	)
{
	STATIC verylong x = 0;

	zintoz(d, &x);
	zadd(a, x, b);
	FREESPACE(x);
}

static void
zaddls(
	verylong a,
	verylong b,
	verylong c
	)
{
/* low level c = a+b, b shorter than a, output can be input */
	register long sa = *a;
	register long sb = *b;
	register long carry = 0;
	register long i;
	register verylong pc;

 /* we know that sa and sb are both >0 or both <0 */
	pc = &c[0];
	if (sa < 0)
	{
		sa = -sa;
		sb = -sb;
	}
	for (i = 1; i <= sb; i++)
	{
		if ((*(++pc) = (*(++a)) + (*(++b)) + carry) < RADIX)
			carry = 0;
		else
		{
			*pc -= RADIX;
			carry = 1;
		}
	}
	for (; i <= sa; i++)
	{
		if ((*(++pc) = (*(++a)) + carry) < RADIX)
			carry = 0;
		else
		{
			*pc -= RADIX;
			carry = 1;
		}
	}
	if (carry)
	{
		c[0] = sa + 1;
		*(++pc) = 1;
	}
	else
		c[0] = sa;
}

void
zsubpos(
	verylong a,
	verylong b,
	verylong *cc
	)
{
	register long sa = a[0];
	register long sb = b[0];
	register long carry = 0;
	register long i;
	register verylong pc;
	verylong c = *cc;

	if (ALLOCATE && !b)
	{
		if (a)
			zcopy(a, cc);
		else
			zzero(cc);
		return;
	}
	if (ALLOCATE && !a)
	{
		zzero(cc);
		return;
	}
	zsetlength(&c, sa/*, "in zsubpos, third argument"*/);
	/* if *cc == a, then nothing will happen */
	/* if *cc == b, then b might point to freed space, so */
	if (b == *cc) b = c;
	*cc = c;
	pc = &c[0];
	for (i = 1; i <= sb; i++)
	{
		if ((*(++pc) = (*(++a)) - (*(++b)) - carry) >= 0)
			carry = 0;
		else
		{
			*pc += RADIX;
			carry = 1;
		};
	}
	for (; i <= sa; i++)
	{
		if ((*(++pc) = (*(++a)) - carry) >= 0)
			carry = 0;
		else
		{
			*pc += RADIX;
			carry = 1;
		};
	}
	i = sa;
	while ((i > 1) && (!(*pc)))
	{
		i--;
		pc--;
	}
	c[0] = i;
}

void
zadd(
	verylong a,
	verylong b,
	verylong *cc
	)
{
	register long sa;
	register long sb;
	register long anegative;
	verylong c = *cc;

	if (ALLOCATE && !a)
	{
		if (b)
			zcopy(b, cc);
		else
			zzero(cc);
		return;
	}
	if (ALLOCATE && !b)
	{
		zcopy(a, cc);
		return;
	}
	if ((anegative = ((sa = a[0]) < 0)) == ((sb = b[0]) < 0))
	{
	/* signs a and b are the same */
		register long i;
		if (anegative)
		{
			sa = -sa;
			sb = -sb;
		}
		zsetlength(&c, (sa > sb ? sa : sb) + 1/*, "in zadd, third argument"*/);
		if (a == *cc) a = c;
		if (b == *cc) b = c;
		*cc = c;
		if (sa == sb)
		{
			register verylong pc;
			pc = &c[0];
			i = 0;
			for (; sa; sa--)
			{
				if ((*(++pc) = (*(++a)) + (*(++b)) + i) < RADIX)
					i = 0;
				else
				{
					*pc -= RADIX;
					i = 1;
				}
			}
			if (i)
			{
				c[0] = sb + 1;
				*(++pc) = 1;
			}
			else
				c[0] = sb;
		}
		else if (sa > sb)
			zaddls(a, b, c);
		else
			zaddls(b, a, c);
		if (anegative)
			c[0] = -c[0];
	/* if anegative, then c cannot be zero */
	}
	else
	{
	/* signs a and b are different */
		verylong old;
		verylong oldc;

		oldc = c;
		if (anegative)
		{
			a[0] = -a[0];
			old = a;
		}
		else
		{
			b[0] = -b[0];
			old = b;
		}
	/* the one that's negative cannot be zero */
		if (!(sa = zcompare(a, b)))
		{
			zzero(&c);
			*cc = c;
			if (anegative)
			{
				if (old != oldc)
					a[0] = -a[0];
			}
			else if (old != oldc)
				b[0] = -b[0];
		}
		else if (sa > 0)
		{
			zsubpos(a, b, &c);
			*cc = c;
			if (anegative)
			{
				c[0] = -c[0];
				if (old != oldc)
					a[0] = -a[0];
			}
			else if (old != oldc)
				b[0] = -b[0];
		}
		else
		{
			zsubpos(b, a, &c);
			*cc = c;
			if (!anegative)
			{
				c[0] = -c[0];
				if (old != oldc)
					b[0] = -b[0];
			}
			else if (old != oldc)
				a[0] = -a[0];
		}
	}
}

void
zsub(
	verylong a,
	verylong b,
	verylong *cc
	)
{
	register long sa;
	register long sb;
	register long anegative;
	verylong c = *cc;

	if (ALLOCATE && !b)
	{
		if (a)
			zcopy(a, cc);
		else
			zzero(cc);
		return;
	}
	if (ALLOCATE && !a)
	{
		zcopy(b, cc);
		znegate(cc);
		return;
	}
	if ((anegative = ((sa = a[0]) < 0)) == ((sb = b[0]) < 0))
	{
	/* signs a and b are the same */
		register long carry = 0;
		register long i;
		register verylong pc;
		register long agrb;

		if (!(agrb = zcompare(a, b)))
		{
			zzero(cc);
			return;
		}
		if ((agrb > 0 && anegative) || (agrb < 0 && !anegative))
		{
			pc = a;
			a = b;
			b = pc;
			sa = *a;
			sb = *b;
		}
		if (anegative)
		{
			sa = -sa;
			sb = -sb;
		}
		zsetlength(&c, sa/*, "in zsub, third argument"*/);
                if (b == *cc) b = c;
		*cc = c;
		pc = &c[0];
		for (i = 1; i <= sb; i++)
		{
			if ((*(++pc) = (*(++a)) - (*(++b)) - carry) >= 0)
				carry = 0;
			else
			{
				*pc += RADIX;
				carry = 1;
			};
		}
		for (; i <= sa; i++)
		{
			if ((*(++pc) = (*(++a)) - carry) >= 0)
				carry = 0;
			else
			{
				*pc += RADIX;
				carry = 1;
			};
		}
		i = sa;
		while ((i > 1) && (!(*pc)))
		{
			i--;
			pc--;
		}
		if (agrb > 0)
			c[0] = i;
		else
			c[0] = -i;
	}
	else
	{
	/* signs a and b are different */
		verylong old;
		verylong oldc;

		oldc = c;
		if (anegative)
		{
			a[0] = -a[0];
			old = a;
		}
		else
		{
			b[0] = -b[0];
			old = b;
		}
	/* the one that's negative cannot be zero */
		zadd(a, b, &c);
		*cc = c;
		if (anegative)
		{
			c[0] = -c[0];
			if (old != oldc)
				a[0] = -a[0];
		}
		else if (old != oldc)
			b[0] = -b[0];
	}
}

void
zsmul(
	verylong a,
	long d,
	verylong *bb
	)
{
	register long sa;
	register long i;
	register verylong pb;
	register long bnegative = 0;
	verylong b = *bb;

	if (( d >= RADIX) || (-d >= RADIX)) {
		STATIC verylong x = 0;
		zintoz(d,&x);
		zmulin(a,&x);
		zcopy(x,bb);
		FREESPACE(x);
		return;
	}

	if (ALLOCATE && !a)
	{
		zzero(bb);
		return;
	}
	if (ALLOCATE && !d)
	{
		zzero(bb);
		return;
	}
	if ((sa = a[0]) < 0)
	{
		sa = (-sa);
		if (d < 0)
			d = (-d);
		else
			bnegative = 1;
	}
	else if ((bnegative = (d < 0)))
		d = (-d);
	zsetlength(&b, sa + 1/*, "in zsmul, third argument"*/);
        if (a == *bb) a = b;
	*bb = b;
	pb = &b[0];
	for (i = sa; i >= 0; i--)
		*pb++ = *a++;
	b[0] = sa;
	sa++;
	*pb = 0;
	zaddmul(d - 1, &b[1], &b[0]);
	while ((sa > 1) && (!(b[sa])))
		sa--;
	b[0] = sa;
	if (bnegative && (b[1] || b[0] != 1))
		b[0] = (-b[0]);
}

static void
kar_mul(
	verylong a,
	verylong b,
	verylong *c,
	long shi
	)
{
	register long al;
	register long hal;
	register long i;
	register long restoreb0 = b[0];
	register verylong pc;
	register long bbig = 1;
	verylong *a0;
	verylong *a1;
	verylong *a2;
	verylong *a3;
	verylong *a4;

	zsetlength(c, (hal = (al = a[0]) + (i = b[0]))/*, "in kar_mul, third argument"*/);
	if ((shi >= (5 * KAR_DEPTH)) || (al < KAR_MUL_CROV) || (i < KAR_MUL_CROV))
	{
		pc = &(*c)[1];
		for (i = hal; i > 0; i--)
			*pc++ = 0;
		pc = &(*c)[1];
		if (al <= *b)
			for (i = al; i; i--)
			{
				zaddmul(*(++a), pc++, b);
			}
		else
			for (i = *b; i; i--)
			{
				zaddmul(*(++b), pc++, a);
			}
		while ((hal > 1) && (!((*c)[hal])))
			hal--;
		(*c)[0] = hal;
		return;
	}
	a0 = &(kar_mem[shi]);
	a1 = &(kar_mem[shi + 1]);
	a2 = &(kar_mem[shi + 2]);
	a3 = &(kar_mem[shi + 3]);
	a4 = &(kar_mem[shi + 4]);
	hal = ((al + 1) >> 1);
	zsetlength(a0, al/*, "in kar_mul, locals\n"*/);
	zsetlength(a1, al/*, ""*/);
	zsetlength(a2, al/*, ""*/);
	zsetlength(a3, al + hal/*, ""*/);
	zsetlength(a4, al + 2/*, ""*/);
	i = hal;
	while ((i > 1) && (!(a[i])))
		i--;
	a[0] = i;
	if (hal >= b[0])
		bbig = 0;
	else
	{
		i = hal;
		while ((i > 1) && (!(b[i])))
			i--;
		b[0] = i;
	}
	for (i = hal + 1; i <= al; i++)
		(*a1)[i - hal] = a[i];
	(*a1)[0] = al - hal;
	if (bbig)
	{
		for (i = hal + 1; i <= restoreb0; i++)
			(*a3)[i - hal] = b[i];
		(*a3)[0] = restoreb0 - hal;
	}
	kar_mul(a, b, a4, shi + 5);
	zadd(a, (*a1), a0);
	a[0] = al;
	if (bbig)
	{
		kar_mul((*a1), (*a3), c, shi + 5);
		zadd(b, (*a3), a2);
		b[0] = restoreb0;
		kar_mul((*a0), (*a2), a3, shi + 5);
	}
	else
		kar_mul((*a0), b, a3, shi + 5);
	zsubpos((*a3), (*a4), a3);
	if (bbig)
		zsubpos((*a3), *c, a3);
	zlshift((*a3), hal * NBITS, a3);
	hal <<= 1;
	if (bbig)
	{
		for (i = (*c)[0]; i; i--)
			(*c)[i + hal] = (*c)[i];
		for (i = hal; i > (*a4)[0]; i--)
			(*c)[i] = 0;
		for (; i; i--)
			(*c)[i] = (*a4)[i];
		(*c)[0] += hal;
	}
	else
	{
		for (i = (*a4)[0]; i >= 0; i--)
			(*c)[i] = (*a4)[i];
	}
	zadd(*c, (*a3), c);
}

void
zmul(
	verylong a,
	verylong b,
	verylong *c
	)
{	/* output not input */
	register long aneg;
	register long bneg;
	verylong olda;
	verylong oldb;

	if (ALLOCATE && (!a || !b))
	{
		zzero(c);
		return;
	}
	if (a == b)
	{
		zsq(a, c);
		return;
	}
	if (!kar_mem_initialized)
	{
		kar_mem_initialized = 1;
		for (aneg = (5 * KAR_DEPTH) - 1; aneg >= 0; aneg--)
			kar_mem[aneg] = (verylong) 0;
	}
	olda = a;
	oldb = b;
	if ((aneg = (*a < 0)))
		a[0] = -a[0];
	if ((bneg = (*b < 0)))
		b[0] = -b[0];
	if (*a > *b)
		kar_mul(a, b, c, (long) 0);
	else
		kar_mul(b, a, c, (long) 0);
	if (aneg != bneg && ((*c)[1] || (*c)[0] != 1))
		(*c)[0] = -(*c)[0];
	if (aneg)
		olda[0] = -olda[0];
	if (bneg)
		oldb[0] = -oldb[0];
}

static void
kar_sq(
	verylong a,
	verylong *c,
	long shi
	)
{
	register long al;
	register long hal;
	register long i;
	register verylong pc;
	verylong *a0;
	verylong *a1;
	verylong *a2;

	zsetlength(c, (i = ((al = a[0]) << 1))/*, "in kar_sq, second argument"*/);
	if ((shi >= (3 * KAR_DEPTH)) || (al < KAR_SQU_CROV))
	{
		register unsigned long uncar;
		long carry = 0;
		pc = &(*c)[1];
		for (; i > 0; i--)
			*pc++ = 0;
		for (hal = 1; hal <= al; hal++)
		{
			i += 2;
			{
				zaddmulsq(al - hal, &((*c)[i]), &(a[hal]));
			}
			uncar = ((*c)[i - 1] << 1) + carry;
			(*c)[i - 1] = uncar & RADIXM;
			uncar = ((*c)[i] << 1) + (uncar >> NBITS);
			{
				zaddmulpsq(&(*c)[i - 1], a[hal], &carry);
			}
			uncar += carry;
			carry = uncar >> NBITS;
			(*c)[i] = uncar & RADIXM;
		}
		while ((i > 1) && (!((*c)[i])))
			i--;
		(*c)[0] = i;
		return;
	}
	a0 = &(kar_mem[shi]);
	a1 = &(kar_mem[shi + 1]);
	a2 = &(kar_mem[shi + 2]);
	hal = ((al + 1) >> 1);
	zsetlength(a0, al + hal + 2/*, "in kar_sq, locals\n"*/);
	zsetlength(a1, al + 2/*, ""*/);
	zsetlength(a2, al/*, ""*/);
	i = hal;
	while ((i > 1) && (!(a[i])))
		i--;
	a[0] = i;
	for (i = hal + 1; i <= al; i++)
		(*a0)[i - hal] = a[i];
	(*a0)[0] = al - hal;
	kar_sq(a, a1, shi + 3);
	zadd(a, (*a0), a2);
	kar_sq((*a0), c, shi + 3);
	a[0] = al;
	kar_sq((*a2), a0, shi + 3);
	zsubpos((*a0), (*a1), a0);
	zsubpos((*a0), *c, a0);
	zlshift((*a0), hal * NBITS, a0);
	hal <<= 1;
	for (i = (*c)[0]; i; i--)
		(*c)[i + hal] = (*c)[i];
	for (i = hal; i > (*a1)[0]; i--)
		(*c)[i] = 0;
	for (; i; i--)
		(*c)[i] = (*a1)[i];
	(*c)[0] += hal;
	zadd(*c, (*a0), c);
}

void
zsq(
	verylong a,
	verylong *c
	)
{	/* output is not input */
	register long aneg;

	if (ALLOCATE && !a)
	{
		zzero(c);
		return;
	}
        if (!kar_mem_initialized)
        {
                kar_mem_initialized = 1;
                for (aneg = (5 * KAR_DEPTH) - 1; aneg >= 0; aneg--)
                        kar_mem[aneg] = (verylong) 0;
        }
	if ((aneg = (*a < 0)))
		a[0] = -a[0];
	kar_sq(a, c, (long) 0);
	if (aneg)
		a[0] = -a[0];
}

void
zmulin(
	verylong b,
	verylong *a
	)
{
	STATIC verylong mem = 0;

	if (ALLOCATE && (!*a || !b))
	{
		zzero(a);
		return;
	}
	zcopy(*a, &mem);
	zmul(mem, b, a);
	FREESPACE(mem);
}

#ifndef ALPHA50
long 
zsdiv(
	verylong a,
	long d,
	verylong *bb
	)
{
	register long sa;
	verylong b = *bb;

#ifndef START
	if (fudge < 0)
		zstart();
#endif
	if (ALLOCATE && !a)
	{
		zzero(bb);
		return (0);
	}
	if (!d)
        {
		zhalt("division by zero in zsdiv");
		return (0);
	}
	if ((sa = a[0]) < 0)
		sa = (-sa);
	zsetlength(&b, sa/*, "in zsdiv, third argument"*/);
	if (a == *bb) a = b;
	*bb = b;
	if ((d >= RADIX) || (d <= -RADIX))
	{
		STATIC verylong zd = 0;
		STATIC verylong zb = 0;
		zintoz(d, &zb);
		zdiv(a, zb, &b, &zd);
		*bb = b;
		sa = ztoint(zd);
		FREE2SPACE(zd,zb);
		return (sa);
	}
	else
	{
		register long den = d;
		register double deninv;
		register long carry = 0;
		register long i;
		long flag = (*a < 0 ? 2 : 0) | (den < 0 ? 1 : 0);

		if (den < 0)
			den = -den;
		deninv = (double)1/den;
		if (a[sa] < den && sa > 1)
			carry = a[sa--];

		for (i = sa; i; i--)
			zdiv21(carry, a[i], den, deninv, b[i]);

		while ((sa > 1) && (!(b[sa])))
			sa--;
		b[0] = sa;
		if (flag)
		{
			if (flag <= 2)
			{
				if (!carry)
					znegate(&b);
				else
				{
					zadd(b, one, &b);
					b[0] = -b[0];
					if (flag == 1)
						carry = carry - den;
					else
						carry = den - carry;
					*bb = b;
				}
			}
			else
				carry = -carry;
		}
		return (carry);
	}
}

#endif


#define correct( q, x1, x2, x3, y1, y2, btopinv) { \
	register long ix1=x1,ix2=x2,ix3=x3,iy1=y1,iy2=y2; \
	long lprodlow=0,lprodhigh=0; \
\
	zaddmulp(&lprodlow, q, iy2, &lprodhigh); \
	if ((ix3 -= lprodlow) < 0) { \
		ix3+= RADIX; \
		ix2--; \
	} \
	lprodlow=0; \
	zaddmulp(&lprodlow, q, iy1, &lprodhigh); \
	if ((ix2 -= lprodlow) < 0) { \
		ix2 += RADIX; \
		--ix1; \
	} \
	if (ix1 < lprodhigh) q--; \
	else if ((ix1 -= lprodhigh)) { \
		q += (long) ((fradix * (fradix * ix1 + ix2))*btopinv ); \
	} \
	else { \
		q += (long) ((fradix * ix2 + ix3)*btopinv ); \
	} \
}


#ifndef ALPHA50
void
zdiv(
	verylong in_a,
	verylong in_b,
	verylong *qqq,
	verylong *rrr
	)

{
	register long i;
	register long qq;
	long sa;
	long sb;
	long sq;
	verylong p;
	verylong pc;
	long sign;
	STATIC verylong a = 0;
	STATIC verylong b = 0;
	STATIC verylong c = 0;
	STATIC verylong d = 0;
	double btopinv;
#ifdef ALPHA
	double btopinv2;
#endif
	double aux;
	verylong q = *qqq;
	verylong r = *rrr;

#ifndef START
	if (fudge < 0)
		zstart();
#endif

	if (ALLOCATE && !in_a)
	{
		zzero(qqq);
		zzero(rrr);
		return;
	}
	if ((!in_b) || (((sb = in_b[0]) == 1) && (!in_b[1])))
	{
		zhalt("division by zero in zdiv");
		return;
	}
	zcopy(in_a,&a);

	zcopy(in_b,&b);

	sign = (*a < 0 ? 2 : 0) | (*b < 0 ? 1 : 0);
	if (*a < 0)
		(*a) = (-(*a));
	sa = (*a);
	if (*b < 0)
		(*b) = (-(*b));
	sb = (*b);
	zsetlength(&c, (i = (sa > sb ? sa : sb) + 1)/*, "in zdiv, locals\n"*/);
	zsetlength(&d, i/*, ""*/);
	zsetlength(&q, i/*, "in zdiv, third argument"*/);
	*qqq = q;
	zsetlength(&r, sb + 1/*, "in zdiv, fourth argument"*/);
	*rrr = r;
	p = &b[sb];
	if ((sb == 1) && (*p < RADIX))
		zintoz(zsdiv(a, *p, &q), &r);
	else
	{
		sq = sa - sb;	/* size of quotient */
		btopinv = (double) (*p) * fradix;
		if (sb > 1)
			btopinv += (*(--p));
		btopinv *= fradix;
#ifdef ALPHA
		btopinv2=btopinv;
		btopinv2 = fudge2/btopinv2;
#else
		if (sb > 2)
			btopinv += (*(p - 1));
#endif
		btopinv = fudge / btopinv;
		p = &a[1];
		pc = &c[0];
		for (i = sa; i > 0; i--)
			*pc++ = *p++;
		p = pc;
		sa = 1 - sb;
		for (i = (-sq); i > 0; i--)
			*p++ = 0;
		*pc = 0;
		d[1] = 0;
		p = &d[sq + 1];
		for (i = sq; i >= 0; i--)
		{

			aux = fradix * (fradix * (*pc) + (*(pc - 1))) + 1.0;
#ifndef ALPHA
			if (i > sa)
				aux += (*(pc - 2));
#endif
#ifdef ALPHA
			qq = (long) (btopinv2 * aux + 1.0);
#else
			qq = (long) (btopinv * aux + 0.5);
		/* dirty, but safe. On most machines (long)(btopinv*aux) */
		/* is correct, or one too big; on some however it becomes */
		/* too small. Could change zstart, but +0.5 and a while */
		/* instead of one single if is safer */
#endif
			if (qq > 0)
			{
#ifdef ALPHA
				if (qq > (1L<<48)) {
					correct(qq,*pc,*(pc-1),(i>sa) ? *(pc-2):0,b[sb],b[sb-1],btopinv);
					if (qq >= RADIX)
						qq = RADIX-1;
				}
#else
				if (qq >= RADIX)
					qq = RADIX-1;
#endif
				zsubmul(qq, &c[i], &b[0]);
				while (*pc < 0)
				{
					qq--;
					{
						zaddmulone(&c[i], &b[0]);
					}
				}
			}
			pc--;
			*p-- = qq;
		}
		sb--;
		while ((sb > 0) && (!(c[sb])))
			sb--;
		sb++;
		r[0] = sb;
		p = &r[1];
		pc = &c[0];
		for (i = sb; i > 0; i--)
			*p++ = *pc++;
		if (sq < 0)
		{
			q[0] = 1;
			q[1] = d[1];
		}
		else
		{
			sq++;
			while ((sq > 1) && (!(d[sq])))
				sq--;
			q[0] = sq;
			p = &q[1];
			pc = &d[1];
			for (i = sq; i > 0; i--)
				*p++ = *pc++;
		}
	}	/* non-trivial case */

	if (sign)
	{
		if (sign <= 2)
		{
			if (!(r[1]) && (r[0] == 1))
				q[0] = -(q[0]);
			else
			{
				zadd(q, one, &q);
				q[0] = -q[0];
				if (sign == 1)
					zsub(r, b, &r);
				else
					zsub(b, r, &r);
			}
		}
		else
			znegate(&r);
	}
	*qqq = q;
	*rrr = r;
	FREE2SPACE(a,b); FREE2SPACE(c,d);
}


void
zmod(
	verylong in_a,
	verylong in_b,
	verylong *rr
	)
{
	register long i;
	register long qq;
	verylong r = *rr;
	STATIC verylong a = 0;
	STATIC verylong b = 0;
	STATIC verylong c = 0;
	long sa;
	long sb;
	long sq;
	verylong p;
	verylong pc;
	long sign;
	double btopinv;
#ifdef ALPHA
        double btopinv2;
#endif
	double aux;

/*printf("in zmod: "); zwrite(in_a); printf(" mod "); zwriteln(in_b); fflush(stdout);
*/
#ifndef START
	if (fudge < 0)
		zstart();
#endif
	if (ALLOCATE && !in_a)
	{
		zzero(rr);
		return;
	}
	if ((!in_b) || (((sb = in_b[0]) == 1) && (!in_b[1])))
	{
		zhalt("division by zero in zmod");
		return;
	}
	zcopy(in_a,&a);
	zcopy(in_b,&b);
	sign = (*a < 0 ? 2 : 0) | (*b < 0 ? 1 : 0);
	if (*a < 0)
		(*a) = (-(*a));
	sa = (*a);
	if (*b < 0)
		(*b) = (-(*b));
	sb = (*b);
	zsetlength(&c, (sa > sb ? sa : sb) + 1/*, "in zmod, local"*/);
	zsetlength(&r, sb + 1/*, "in zmod, third argument"*/);
	*rr = r;
	p = &b[sb];
	if ((sb == 1) && (*p < RADIX))
		zintoz(zsdiv(a, *p, &c), &r);
	else
	{
		sq = sa - sb;
		btopinv = (double) (*p) * fradix;
		if (sb > 1)
			btopinv += (*(--p));
		btopinv *= fradix;
#ifdef ALPHA
                btopinv2=btopinv;
		btopinv2 = fudge2/btopinv2;
#else
                if (sb > 2)
                        btopinv += (*(p - 1));
#endif
		btopinv = fudge / btopinv;
		p = &a[1];
		pc = &c[0];
		for (i = sa; i > 0; i--)
			*pc++ = *p++;
		p = pc;
		sa = 1 - sb;
		for (i = (-sq); i > 0; i--)
			*p++ = 0;
		*pc = 0;
		for (i = sq; i >= 0; i--)
		{
			aux = fradix * (fradix * (*pc) + (*(pc - 1))) + 1.0;
#ifdef ALPHA
			qq = (long) (btopinv2 * aux + 1.0);
#else
			if (i > sa)
				aux += (*(pc - 2));
			qq = (long) (btopinv * aux + 0.5);

		/* see comment in zdiv */
#endif
			if (qq > 0)
			{
#ifdef ALPHA
				if (qq > (1L<<48)) {
					correct(qq,*pc,*(pc-1),(i>sa) ? *(pc-2):0,b[sb],b[sb-1],btopinv);
					if (qq >= RADIX)
						qq = RADIX-1;
				}
#else
				if (qq >= RADIX)
					qq = RADIX-1;
#endif
				zsubmul(qq, &c[i], &b[0]);
				while (*pc < 0)
				{
					{
						zaddmulone(&c[i], &b[0]);
					}
				}
			}
			pc--;
		}	/* loop on i */
		sb--;
		while ((sb > 0) && (!(c[sb])))
			sb--;
		sb++;
		r[0] = sb;
		p = &r[1];
		pc = &c[0];
		for (i = sb; i > 0; i--)
			*p++ = *pc++;
	}	/* non-trivial case */

	if (sign)
	{
		if ((sign <= 2) && (!((r[0] == 1) && !(r[1]))))
		{
			if (sign == 1)
				zsub(r, b, &r);
			else
				zsub(b, r, &r);
		}
		else
			znegate(&r);	/* negating r==0 doesn't hurt */
	}
	*rr = r;
	FREE2SPACE(a,b); FREESPACE(c);

}
#endif


void
zmulmod(
	verylong a,
	verylong b,
	verylong n,
	verylong *c
	)
{
	STATIC verylong mem = 0;

	if (ALLOCATE && !n)
	{
		zhalt("modulus zero in zmulmod");
		return;
	}
	if (ALLOCATE && (!a || !b))
	{
		zzero(c);
		return;
	}
	zsetlength(&mem, a[0] + b[0] + 2/*, "in zmulmod, local"*/);
	zmul(a, b, &mem);
	zmod(mem, n, c);
	FREESPACE(mem);
}

void
zsqmod(
	verylong a,
	verylong n,
	verylong *c
	)
{
	STATIC verylong mem = 0;

	if (ALLOCATE && !n)
	{
		zhalt("modulus zero in zsqmod");
		return;
	}
	if (ALLOCATE && !a)
	{
		zzero(c);
		return;
	}
	zsetlength(&mem, 2 * a[0] + 2/*, "in zsqmod, local"*/);
	zsq(a, &mem);
	zmod(mem, n, c);
	FREESPACE(mem);
}

static long 
zxxeucl(
	verylong ain,
	verylong nin,
	verylong *invv,
	verylong *uu
	)
{
	STATIC verylong a = 0;
	STATIC verylong n = 0;
	STATIC verylong q = 0;
	STATIC verylong w = 0;
	STATIC verylong x = 0;
	STATIC verylong y = 0;
	STATIC verylong z = 0;
	verylong inv = *invv;
	verylong u = *uu;
	long diff;
	long ilo;
	long sa;
	long sn;
	long temp;
	long e;
	long fast;
	long parity;
	long gotthem;
	verylong pin;
	verylong p;
	long i;
	long try11;
	long try12;
	long try21;
	long try22;
	long got11;
	long got12;
	long got21;
	long got22;
	double hi;
	double lo;
	double dt;
	double fhi;
	double flo;
	double num;
	double den;
	double dirt;

#ifndef START
	if (fudge < 0)
		zstart();
#endif
	zsetlength(&a, (e = (ain[-1] > nin[-1] ? ain[-1] : nin[-1]))/*, 
	           "in zxxeucl, locals\n"*/);
	zsetlength(&n, e/*, ""*/);
	zsetlength(&q, e/*, ""*/);
	zsetlength(&w, e/*, ""*/);
	zsetlength(&x, e/*, ""*/);
	zsetlength(&y, e/*, ""*/);
	zsetlength(&z, e/*, ""*/);
	pin = &ain[0];
	p = &a[0];
	for (i = (*pin); i >= 0; i--)
		*p++ = *pin++;
	zsetlength(&inv, e/*, "in zxxeucl, third argument"*/);
	*invv = inv;
	zsetlength(&u, e/*, "in zxxeucl, fourth argument"*/);
	*uu = u;
	fhi = 1.0 + epsilon;
	flo = 1.0 - epsilon;
	pin = &nin[0];
	p = &n[0];
	for (i = (*pin); i >= 0; i--)
		*p++ = *pin++;
	inv[0] = 1;
	inv[1] = 1;
	w[0] = 1;
	w[1] = 0;
	while (n[0] > 1 || n[1] > 0)
	{
		gotthem = 0;
		sa = a[0];
		sn = n[0];
		diff = sa - sn;
		if (!diff || diff == 1)
		{
			sa = a[0];
			p = &a[sa];
			num = (double) (*p) * fradix;
			if (sa > 1)
				num += (*(--p));
			sn = n[0];
			p = &n[sn];
			den = (double) (*p) * fradix;
			if (sn > 1)
				den += (*(--p));
 			hi = fhi * (num + 1.0) / den;
			lo = flo * num / (den + 1.0);
			if (diff > 0)
			{
				hi *= fradix;
				lo *= fradix;
			}
			try11 = 1;
			try12 = 0;
			try21 = 0;
			try22 = 1;
			parity = 1;
			fast = 1;
			while (fast > 0)
			{
				parity = 1 - parity;
				if (hi >= fradix)
					fast = 0;
				else
				{
					ilo = (long)lo;
					dirt = hi - ilo;
					if (!dirt || !ilo || ilo < (long)hi)
						fast = 0;
					else
					{
						dt = lo;
						lo = flo / dirt;
						if (dt > ilo)
							hi = fhi / (dt - ilo);
						else
							hi = fradix;
						temp = try11;
						try11 = try21;
						if ((RADIX - temp) / ilo < try21)
							fast = 0;
						else
							try21 = temp + ilo * try21;
						temp = try12;
						try12 = try22;
						if ((RADIX - temp) / ilo < try22)
							fast = 0;
						else
							try22 = temp + ilo * try22;
						if ((fast > 0) && (parity > 0))
						{
							gotthem = 1;
							got11 = try11;
							got12 = try12;
							got21 = try21;
							got22 = try22;
						}
					}
				}
			}
		}
		if (gotthem)
		{
			zsmul(inv, got11, &x);
			zsmul(w, got12, &y);
			zsmul(inv, got21, &z);
			zsmul(w, got22, &w);
			zadd(x, y, &inv);
			zadd(z, w, &w);
			zsmul(a, got11, &x);
			zsmul(n, got12, &y);
			zsmul(a, got21, &z);
			zsmul(n, got22, &n);
			zsub(x, y, &a);
			zsub(n, z, &n);
		}
		else
		{
			zdiv(a, n, &q, &a);
			zmul(q, w, &x);
			zadd(inv, x, &inv);
			if (a[0] > 1 || a[1] > 0)
			{
				zdiv(n, a, &q, &n);
				zmul(q, inv, &x);
				zadd(w, x, &w);
			}
			else
			{
				p = &a[0];
				pin = &n[0];
				for (i = (*pin); i >= 0; i--)
					*p++ = *pin++;
				n[0] = 1;
				n[1] = 0;
				zsub(nin, w, &inv);
			}
		}
	}
	if ((a[0] == 1) && (a[1] == 1))
		e = 0;
	else
		e = 1;
	p = &u[0];
	pin = &a[0];
	for (i = (*pin); i >= 0; i--)
		*p++ = *pin++;
	*invv = inv;
	*uu = u;
	FREE2SPACE(a,n); FREE2SPACE(q,w); FREE2SPACE(x,y); FREESPACE(z);
	return (e);
}

long 
zinv(
	verylong ain,
	verylong nin,
	verylong *invv
	)
{
	STATIC verylong u = 0;

	if ((!nin) || (!ain) || (ain[0] < 0) || (nin[0] < 0))
	{
		zhalt("zero or negative argument(s) in zinv");
		return (0);
	}
	if (!(zxxeucl(ain, nin, invv, &u))) {
		FREESPACE(u);
		return (0);
	}
	zcopy(u, invv);
	FREESPACE(u);
	return (1);
}

void
zexpmod(
	verylong a,
	verylong e,
	verylong n,
	verylong *bb
	)
{
	register long i;
	register long j;
	register long k = 0;
	verylong b = *bb;
	STATIC verylong loca = 0;

	if (ALLOCATE && !n)
	{
		zhalt("modulus zero in zexpmod");
		return;
	}
	if (ALLOCATE && !e)
	{
		zone(bb);
		return;
	}
	if (ALLOCATE && !a)
	{
		zzero(bb);
		return;
	}
	if ((i = e[0]) < 0)
		i = (-i);
	if ((i == 1) && (!(e[1])))
		zone(&b);
	else
	{
		zmod(a, n, &loca);
		zcopy(loca, &b);
		for (; i; i--)
		{
			j = e[i];
			if (!k)
			{
				k = 1;
				while ((k << 1) <= j)
					k <<= 1;
			}
			while (k >>= 1)
			{
				zsqmod(b, n, &b);
				if (j & k)
					zmulmod(loca, b, n, &b);
			}
			k = RADIX;
		}
	}
	if (e[0] < 0 && zinv(b, n, &b))
		zhalt("undefined quotient in zexpmod");
	*bb = b;
	FREESPACE(loca);
}

void
z2mul(
	verylong n,
	verylong *rres
	)
{
	register long sn;
	register long i;
	register long carry = 0;
	verylong res = *rres;

	if (ALLOCATE && !n)
	{
		zzero(rres);
		return;
	}
	if ((!n[1]) && (n[0] == 1))
	{
		zzero(rres);
		return;
	}
	if ((sn = n[0]) < 0)
		sn = -sn;
	zsetlength(&res, sn + 1/*, "in z2mul, second argument"*/);
	if (n == *rres) n = res;
	*rres = res;
	for (i = 1; i <= sn; i++)
	{
		if ((res[i] = (n[i] << 1) + carry) >= RADIX)
		{
			res[i] -= RADIX;
			carry = 1;
		}
		else
			carry = 0;
	}
	if (carry)
		res[++sn] = 1;
	if (n[0] < 0)
		res[0] = -sn;
	else
		res[0] = sn;
}

long 
z2div(
	verylong n,
	verylong *rres
	)
{
	register long sn;
	register long i;
	register long result;
	verylong res = *rres;

	if ((!n) || ((!n[1]) && (n[0] == 1)))
	{
		zzero(rres);
		return (0);
	}
	if ((sn = n[0]) < 0)
		sn = -sn;
	zsetlength(&res, sn/*, "in z2div, second argument"*/);
	if (n == *rres) n = res;
	*rres = res;
	result = n[1] & 1;
	for (i = 1; i < sn; i++)
	{
		res[i] = (n[i] >> 1);
		if (n[i + 1] & 1)
			res[i] += (RADIX >> 1);
	}
	if ((res[sn] = (n[sn] >> 1)))
		res[0] = n[0];
	else if (sn == 1)
	{
		res[0] = 1;
	}
	else if (n[0] < 0)
		res[0] = -sn + 1;
	else
		res[0] = sn - 1;
	return (result);
}

long 
z2mod(
	verylong n
	)
{
	if (!n)
		return (0);
	return (n[1] & 1);
}

void
zlshift(
	verylong n,
	long k,
	verylong *rres
	)
{
	register long big;
	register long smallval;
	register long sn;
	register long i;
	register long cosmall;
	verylong res = *rres;

	if (ALLOCATE && !n)
	{
		zzero(rres);
		return;
	}
	if ((!n[1]) && (n[0] == 1))
	{
		zzero(rres);
		return;
	}
	if (!k)
	{
		if (n != *rres)
			zcopy(n, rres);
		return;
	}
	if (k < 0)
	{
		zrshift(n, -k, rres);
		return;
	}
	if (k == 1)
	{
		z2mul(n, rres);
		return;
	}
	if ((sn = n[0]) < 0)
		sn = -sn;
	i = sn + (big = k / NBITS);
	if ((smallval = k - big * NBITS))
	{
		zsetlength(&res, i + 1/*, "in zlshift, third argument"*/);
		if (n == *rres) n = res;
		*rres = res;
		res[i + 1] = n[sn] >> (cosmall = NBITS - smallval);
		for (i = sn; i > 1; i--)
			res[i + big] = ((n[i] << smallval) & RADIXM) + (n[i - 1] >> cosmall);
		res[big + 1] = (n[1] << smallval) & RADIXM;
		for (i = big; i; i--)
			res[i] = 0;
		if (res[sn + big + 1])
			big++;
	}
	else
	{
		zsetlength(&res, i/*, "in zlshift, third argument"*/);
		if (n == *rres) n = res;
		*rres = res;
		for (i = sn; i; i--)
			res[i + big] = n[i];
		for (i = big; i; i--)
			res[i] = 0;
	}
	if (n[0] > 0)
		res[0] = n[0] + big;
	else
		res[0] = n[0] - big;
}

void
zrshift(
	verylong n,
	long k,
	verylong *rres
	)
{
	register long big;
	register long smallval;
	register long sn;
	register long i;
	register long cosmall;
	verylong res = *rres;

	if (ALLOCATE && !n)
	{
		zzero(rres);
		return;
	}
	if ((!n[1]) && (n[0] == 1))
	{
		zzero(rres);
		return;
	}
	if (!k)
	{
		if (n != *rres)
			zcopy(n, rres);
		return;
	}
	if (k < 0)
	{
		zlshift(n, -k, rres);
		return;
	}
	if (k == 1)
	{
		z2div(n, rres);
		return;
	}
	big = k / NBITS;
	smallval = k - big * NBITS;
	if ((sn = n[0]) < 0)
		sn = -sn;
	if ((big >= sn) || ((big == sn - 1) && (!(n[sn] >> smallval))))
	{
		zzero(rres);
		return;
	}
	sn -= big;
	zsetlength(&res, sn/*, "in zrshift, third argument"*/);
	if (n == *rres) n = res;
	*rres = res;
	if (smallval)
	{
		cosmall = NBITS - smallval;
		for (i = 1; i < sn; i++)
			res[i] = (n[i + big] >> smallval) +
				((n[i + big + 1] << cosmall) & RADIXM);
		if (!(res[sn] = (n[sn + big] >> smallval)))
			sn--;
	}
	else
		for (i = 1; i <= sn; i++)
			res[i] = n[i + big];
	if (n[0] > 0)
		res[0] = sn;
	else
		res[0] = -sn;
}

long 
value(
	char c
	)
{
	switch (c)
	{
		case '0':
			return (0);
		case '1':
			return (1);
		case '2':
			return (2);
		case '3':
			return (3);
		case '4':
			return (4);
		case '5':
			return (5);
		case '6':
			return (6);
		case '7':
			return (7);
		case '8':
			return (8);
		case '9':
			return (9);
		case 'A':
		case 'a':
			return (10);
		case 'B':
		case 'b':
			return (11);
		case 'C':
		case 'c':
			return (12);
		case 'D':
		case 'd':
			return (13);
		case 'E':
		case 'e':
			return (14);
		case 'F':
		case 'f':
			return (15);
		default:
			return(0);
	}
}

void 
zhsread(
	char *str,
	verylong *aa
	)
{
	// static char *inmem = 0;
	char *inmem;
	char *in;
	register long d = 0;
	register long anegative = 0;
	verylong a = *aa;
	extern double log();

	if (log16rad < 0)
		log16rad = log((double)RADIX) / log((double)16);

	//if (!inmem)
		inmem = (char *)calloc((size_t)IN_LINE, sizeof(char));
	if (sscanf(str, "%s", inmem) == EOF)
    {
        free (inmem);
		return;
    }
	if ((inmem[0] == '-') || (inmem[0] == '_'))
	{
		anegative = 1;
		in = (inmem + 1);
	}
	else
		in = inmem;
	zsetlength(&a, (long) (strlen(in) / log16rad)/*, "in zhsread, second argument"*/);
	a[0] = 1;
	a[1] = 0;
	while (in[d] != '\0')
	{
		if (in[d] == IN_LINE_BREAK)
		{
			if (sscanf(str, "%s", in) == EOF)
            {
                free (inmem);
				return;
            }
			d = 0;
			if (strlen(in) > log16rad * (a[-1] - a[0]))
			{
				zsetlength(&a, (long) (a[0] + 3 + strlen(in) / log16rad)/*, 
				           "in zhsread, second argument"*/);
			}
		}
		else
		{
			zlshift(a, (long) 4, &a);
			zsadd(a,value(in[d++]),&a);
		}
	}
	if (anegative)
		znegate(&a);
	*aa = a;

    free (inmem);
}

#ifdef ALPHA50
#define lower_radix(in,out) { \
register long i=1,j=0,rb=0,nrb=0,bits; \
register verylong pin=in; \
long size; \
 \
if ((size= *pin)==2) { \
	pin++; \
	out[1] = *pin & ALPHA50RADIXM; \
	out[2] = *pin >> ALPHA50NBITS; \
	pin++; \
	out[2] |= (*pin & ((1L << (2*ALPHA50NBITS-NBITS))-1L)) << (NBITS - \
	 ALPHA50NBITS); \
	bits = *pin >> (2*ALPHA50NBITS-NBITS); \
	if (bits) { \
	        out[3] = bits; \
	        out[0]=3; \
	        } \
	else \
	        out[0]=2; \
	} \
else { \
	for (j=0; i <= size ; j++) { \
	        if (nrb >= ALPHA50NBITS) { \
	                bits = rb & ALPHA50RADIXM; \
	                rb >>= ALPHA50NBITS; \
	                nrb -= ALPHA50NBITS; \
	                } \
	        else { \
	                bits = rb | ((((1L << (ALPHA50NBITS - nrb)) -1L) & \
	                 pin[i]) << nrb); \
	                rb = pin[i] >> (ALPHA50NBITS - nrb); \
	                nrb += (NBITS-ALPHA50NBITS); \
	                i++; \
	                } \
	        out[j+1]=bits; \
	        } \
	if (rb) { \
	        j++; \
	        if (rb < (1L<<ALPHA50NBITS)) \
	                out[j]= rb; \
	        else { \
	                out[j]= rb & (( 1L << ALPHA50NBITS) -1); \
	                j++; \
	                out[j] = rb>>ALPHA50NBITS; \
	                } \
	        } \
	out[0]=j; \
	} \
}


/* convert from ALPHA50NBITS to NBITS radix */
/* this can probably be made faster */

#define raise_radix(size,in,out) { \
register long i,j,nrb,bits; \
 \
nrb=ALPHA50NBITS; \
for (i=0,j=1; i<size; j++) { \
	bits = (in[i++] >> (ALPHA50NBITS-nrb)); \
	if (i <size ) { \
	        if (nrb + ALPHA50NBITS >= NBITS) { \
	                bits = bits | ((in[i] & ( (1L<< (NBITS-nrb)) - 1)) << \
	                 nrb); \
	                nrb -= (NBITS-ALPHA50NBITS); \
	                if (!nrb) { i++; nrb=ALPHA50NBITS; } \
	                } \
	        else { \
	                bits = bits | (in[i++] << nrb); \
	                if (i< size) \
	                        bits = bits | ((in[i] & ((1L <<( (NBITS- \
	                         ALPHA50NBITS)-nrb)) -1)) << (ALPHA50NBITS \
	                         + nrb)); \
	                nrb = 2*ALPHA50NBITS-NBITS+nrb; \
	                if (!nrb) { i++; nrb=ALPHA50NBITS; } \
	                } \
	        } \
	out[j]=bits; \
	} \
if (j>2 && !out[j-1]) j--; \
out[0]=j-1; \
}


#define raise_radix2(size,in,out) { \
register long i,j,nrb,bits; \
 \
nrb=ALPHA50NBITS; \
for (i=1,j=1; i<=size; j++) { \
	bits = (in[i++] >> (ALPHA50NBITS-nrb)); \
	if (i <=size ) { \
	        if (nrb + ALPHA50NBITS >= NBITS) { \
	                bits = bits | ((in[i] & ( (1L<< (NBITS-nrb)) - 1)) << nrb); \
	                nrb -= (NBITS-ALPHA50NBITS); \
	                if (!nrb) { i++; nrb=ALPHA50NBITS; } \
	                } \
	        else { \
	                bits = bits | (in[i++] << nrb); \
	                if (i<=size) \
	                        bits = bits | ((in[i] & ((1L <<( (NBITS- \
	                         ALPHA50NBITS)-nrb)) -1)) << (ALPHA50NBITS+nrb)); \
	                nrb = 2*ALPHA50NBITS-NBITS+nrb; \
	                if (!nrb) { i++; nrb=ALPHA50NBITS; } \
	                } \
	        } \
	out[j]=bits; \
	} \
if (j>2 && !out[j-1]) j--; \
out[0]= (j > 1) ? j-1:1 ; \
}



/* 50 bit division for the alpha */

void
zdiv(
	verylong in_a,
	verylong in_b,
	verylong *qqq,
	verylong *rrr
	)

{
	register long i;
	register long qq;
	long sa;
	long sb;
	long sq;
	verylong p;
	verylong pc;
	long sign;
	STATIC verylong a = 0;
	STATIC verylong b = 0;
	STATIC verylong c = 0;
	STATIC verylong d = 0;
	double btopinv;
	double aux;
	verylong q = *qqq;
	verylong r = *rrr;

/*printf("zdiv: "); zwrite(in_a); printf(" div "); zwriteln(in_b);*/
#ifndef START
	if (fudge < 0)
	        zstart();
#endif
	if (ALLOCATE && !in_a)
	{
	        zzero(qqq);
	        zzero(rrr);
	        return;
	}
	if ((!in_b) || (((sb = in_b[0]) == 1) && (!in_b[1])))
	{
	        zhalt("division by zero in zdiv");
	        return;
	}

	sign = (*in_a < 0 ? 2 : 0) | (*in_b < 0 ? 1 : 0);
	if (*in_a < 0)
	        (*in_a) = (-(*in_a));
	if (*in_b < 0)
	        (*in_b) = (-(*in_b));

	if ((*in_b == 1) && (in_b[1] < ALPHA50RADIX)) {
	        zintoz(zsdiv(in_a, in_b[1], qqq), rrr);
	        goto done;
	        }

	zsetlength(&a,(long) (NBITS/(double)ALPHA50NBITS*(double)in_a[0]+1.0)/*,
	        ""*/);
	lower_radix(in_a,a);    /* test this for large integers */

	/* convert b */
	zsetlength(&b,(long) (NBITS/(double)ALPHA50NBITS*(double)in_b[0]+1.0)/*,
	        ""*/);
	lower_radix(in_b,b);

	sa = (*a);
	sb = (*b);
	zsetlength(&c, (i = (sa > sb ? sa : sb) + 1)/*, "in zdiv, locals\n"*/);
	zsetlength(&d, i/*, ""*/);
	zsetlength(&q, i/*, "in zdiv, third argument"*/);
	*qqq = q;
	zsetlength(&r, sb + 1/*, "in zdiv, fourth argument"*/);
	*rrr = r;
	p = &b[sb];

	sq = sa - sb;
	btopinv = (double) (*p) * alpha50fradix;
	if (sb > 1)
	        btopinv += (*(--p));
	btopinv *= (double) alpha50fradix;
	if (sb > 2)
	        btopinv += (*(p - 1));
	btopinv = alpha50fudge / btopinv;
	p = &a[1];
	pc = &c[0];
	for (i = sa; i > 0; i--)
	        *pc++ = *p++;
	p = pc;
	sa = 1 - sb;
	for (i = (-sq); i > 0; i--)
	        *p++ = 0;
	*pc = 0;
	d[1] = 0;
	p = &d[sq + 1];
	for (i = sq; i >= 0; i--)
	{
	        aux = alpha50fradix * (alpha50fradix * (*pc) +
	            (*(pc - 1))) + 1.0;
	        if (i > sa)
	                aux += (*(pc - 2));
	        qq = (long) (btopinv * aux + 0.5);
	/* dirty, but safe. On most machines (long)(btopinv*aux) */
	/* is correct, or one too big; on some however it becomes */
	/* too small. Could change zstart, but +0.5 and a while */
	/* instead of one single if is safer */
	        if (qq >= ALPHA50RADIX)
	                qq = ALPHA50RADIX-1;
	        if (qq > 0)
	        {
	                zalphasubmul(qq, &c[i], &b[0]);
	                while (*pc < 0)
	                {
	                        qq--;
	                        {
	                                zalphaaddmulone(&c[i], &b[0]);
	                        }
	                }
	        }
	        pc--;
	        *p-- = qq;
	}
	sb--;
	while ((sb > 0) && (!(c[sb])))
	        sb--;
	sb++;
	raise_radix(sb,c,r);

	if (sq < 0)
	{
	        q[0] = 1;
	        q[1] = d[1];
	}
	else
	{
	        sq++;
	        while ((sq > 1) && (!(d[sq])))
	                sq--;
	        raise_radix2(sq,d,q);
	}

	if (sign)
	{
	        if (sign <= 2)
	        {
	                if (!(r[1]) && (r[0] == 1))
	                        q[0] = -(q[0]);
	                else
	                {
	                        zadd(q, one, &q);
	                        q[0] = -q[0];
	                        if (sign == 1)
	                                zsub(r, b, &r);
	                        else
	                                zsub(b, r, &r);
	                }
	        }
	        else
	                znegate(&r);
	}
	*qqq = q;
	*rrr = r;
done:	;

	FREE2SPACE(a,b); FREE2SPACE(c,d);
}


void
zmod(
	verylong in_a,
	verylong in_b,
	verylong *rr
	)
{
	register long i;
	register long qq;
	verylong r = *rr;
	STATIC verylong a = 0;
	STATIC verylong b = 0;
	STATIC verylong c = 0;
	long sa;
	long sb;
	long sq;
	verylong p;
	verylong pc;
	long sign;
	double btopinv;
	double aux;

#ifndef START
	if (fudge < 0)
	        zstart();
#endif
	if (ALLOCATE && !in_a)
	{
	        zzero(rr);
	        return;
	}
	if ((!in_b) || (((sb = in_b[0]) == 1) && (!in_b[1])))
	{
	        zhalt("division by zero in zmod");
	        return;
	}

	sign = (*in_a < 0 ? 2 : 0) | (*in_b < 0 ? 1 : 0);
	if (*in_a < 0)
	        (*in_a) = (-(*in_a));
	if (*in_b < 0)
	        (*in_b) = (-(*in_b));

	if ((*in_b == 1) && (in_b[1] < ALPHA50RADIX)) {
	        zintoz(zsdiv(in_a, in_b[1], &c), rr);
	        goto done;
	        }

	  zsetlength(&a,(long) (NBITS/(double)ALPHA50NBITS*(double)in_a[0]+1.0)
            /*, ""*/);
	  lower_radix(in_a,a);    /* test this for large integers */

	  zsetlength(&b,(long) (NBITS/(double)ALPHA50NBITS*(double)in_b[0]+1.0)
            /*, ""*/);
	  lower_radix(in_b,b);

	sa = (*a);
	sb = (*b);
	zsetlength(&c, (sa > sb ? sa : sb) + 1/*, "in zmod, local"*/);
	zsetlength(&r, sb + 1/*, "in zmod, third argument"*/);
	*rr = r;
	p = &b[sb];

	sq = sa - sb;
	btopinv = (double) (*p) * alpha50fradix;
	if (sb > 1)
	        btopinv += (*(--p));
	btopinv *= alpha50fradix;
	if (sb > 2)
	        btopinv += (*(p - 1));
	btopinv = alpha50fradix / btopinv;
	p = &a[1];
	pc = &c[0];
	for (i = sa; i > 0; i--)
	        *pc++ = *p++;
	p = pc;
	sa = 1 - sb;
	for (i = (-sq); i > 0; i--)
	        *p++ = 0;
	*pc = 0;
	for (i = sq; i >= 0; i--)
	{
	        aux = alpha50fradix * (alpha50fradix * (*pc) +
	         (*(pc - 1))) + 1.0;
	        if (i > sa)
	                aux += (*(pc - 2));
	        qq = (long) (btopinv * aux + 0.5);
	/* see comment in zdiv */
	        if (qq >= ALPHA50RADIX)
	                qq = ALPHA50RADIX-1;
	        if (qq > 0)
	        {
	                zalphasubmul(qq, &c[i], &b[0]);
	                while (*pc < 0)
	                {
	                        {
	                                zalphaaddmulone(&c[i], &b[0]);
	                        }
	                }
	        }
	        pc--;
	}       /* loop on i */
	sb--;
	while ((sb > 0) && (!(c[sb])))
	        sb--;
	sb++;
	raise_radix(sb,c,r);

	if (sign)
	{
	        if ((sign <= 2) && (!((r[0] == 1) && !(r[1]))))
	        {
	                if (sign == 1)
	                        zsub(r, b, &r);
	                else
	                        zsub(b, r, &r);
	        }
	        else
	                znegate(&r);    /* negating r==0 doesn't hurt */
	}
	*rr = r;
done:	;
	FREE2SPACE(a,b); FREESPACE(c);
}


static void
zalphasubmul(
	  long r,
	  verylong a,
	  verylong b
	  )
{
	  register long rd = ALPHA50RADIX - r;
	  register long i;
	  long carry = ALPHA50RADIX;

	  for (i = (*b++); i > 0; i--)
	  {
	          zaldivaddmulp(a, *b, rd, &carry);
	          a++;
	          carry += ALPHA50RADIXM - (*b++);
	  }
	  *a += carry - ALPHA50RADIX;    /* unnormalized */
}




/* alpha zsdiv */

long 
zsdiv(
	verylong in_a,
	long d,
	verylong *bb
	)
{
	register long sa;
	verylong b = *bb;
	STATIC verylong a=0;
	long flag;

#ifndef START
	if (fudge < 0)
	        zstart();
#endif
	if (ALLOCATE && !in_a)
	{
	        zzero(bb);
	        return (0);
	}
	if (!d)
	{
	        zhalt("division by zero in zsdiv");
	        return (0);
	}

	if ((d >= ALPHA50RADIX) || (d <= -ALPHA50RADIX))
	{
	        STATIC verylong zd = 0;
	        STATIC verylong zb = 0;
		long keepzd;

	        zintoz(d, &zb);
	        zdiv(in_a, zb, bb, &zd);
		keepzd = ztoint(zd);
		FREE2SPACE(zd,zb);
	        return (keepzd);
	}

	if (*in_a < 0) {
	        *in_a = - (*in_a);
	        flag = 2L |  (d < 0 ? 1L : 0L);
	        }
	else
	        flag = (d < 0 ? 1L : 0L);

	zsetlength(&a,(long) (NBITS/(double)ALPHA50NBITS*(double) *in_a+1.0)/*, ""*/);
	lower_radix(in_a,a);
	sa = a[0];

	zsetlength(&b, sa/*, "in zsdiv, third argument"*/);
	*bb = b;
	{
	        register long den = d;
	        register double deninv;
	        register long carry = 0;
	        register long i;

	        if (den < 0)
	                den = -den;
	        deninv = (double)1.0/den;
	        if (a[sa] < den && sa > 1)
	                carry = a[sa--];
	        for (i = sa; i; i--)
	        {
	                zdiv21(carry, a[i], den, deninv, b[i]);
	        }
	        while ((sa > 1) && (!(b[sa])))
	                sa--;
	        b[0] = sa;

	        raise_radix2(sa,b,(*bb));

	        if (flag)
	        {
	                if (flag <= 2)
	                {
	                        if (!carry)
	                                znegate(bb);
	                        else
	                        {
	                                zadd(*bb, one, bb);
	                                *bb[0] = -(*bb[0]);
	                                if (flag == 1)
	                                        carry = carry - den;
	                                else
	                                        carry = den - carry;
	                        }
	                }
	                else
	                        carry = -carry;
	        }
		FREESPACE(a);
	        return (carry);
	}
	FREESPACE(a);
}


#endif

#endif

